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ABSTRACT 

Context. When trying to derive the star cluster physical parameters of the M33 galaxy using broad-band unresolved ground-based 
photometry, previous studies mainly made use of simple stellar population models, shown in the recent years to be oversimplified. 
Aims. In this study, we aim to derive the star cluster physical parameters (age, mass, and extinction; metallicity is assumed to be 
LMC-like for clusters with age below 1 Gyr and left free for older clusters) of this galaxy using models that take stochastic dispersion 
of cluster integrated colors into account. 

Methods. We use three recently published M33 catalogs of cluster optical broad-band photometry in standard UBVRI and in 
CFHT/MegaCam u*g'r'i'z' photometric systems. We also use near-infrared JHK photometry that we derive from deep 2MASS 
images. We derive the cluster parameters using a method that takes into account the stochasticity problem, presented in previous 
papers of this series. 

Results. The derived differential age distribution of the M33 cluster population is composed of a two-slope profile indicating that 
the number of clusters decreases when age gets older. The first slope is interpreted as the evolutionary fading phase of the cluster 
magnitudes, and the second slope as the cluster disruption. The threshold between these two phases occurs at ~300Myrs, comparable 
to what is observed in the M31 galaxy. We also model by use of artificial clusters the ability of the cluster physical parameter derivation 
method to correctly derive the two-slope profile for different photometric systems tested. 
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1. Introduction 


There is a current need for an accurate catalog for the star clus¬ 
ter system of the Triangulum galaxy, or Messier 33 (M33), as it 
could be used as a constraint for the derivation of the galaxy star 
formation history. Several other reasons encourage the study of 
this particular star cluster system. The nearly face-on inclination 
(i = 56 degrees, [Regan & Vogel|1994) of M33 reduces extinction 


effects for the majority of its c 


uster population, situated in the 


disk. Also, M33 is the only close late-type spiral galaxy, situated 
at a distance of 867 kpc ( Ga lleti et al.||2004| distance modulus 
of (m - M )o = 24.69), making its star cluster system accessi¬ 
ble to ground-based telescopes for integrated photometric and 
spectroscopic studies and to the Hubble Space Telescope (HST) 
for resolved measurements. While other star cluster systems of 
Local Group galaxies have received considerable attention, as in 
the case of M31 and the Magellanic Clouds, the M33 star clus¬ 
ter system has not been studied as much. Therefore an extended 
knowledge of its star cluster system would improve the under¬ 
standing of the relationship between star clusters and their host 
galaxies. 

The M33 star cluster system has nevertheless been studied 


1960 

Melnick & D’Odorico 

1978; Christian & Schommer 1982, 

1988 

Mochejskaet al. 1998 

contributed to the building of a cat- 


alog of star clusters in M3 3 using ground-based unresolved pho¬ 
tometry in optical passbands. Chan dar et al.| 


used the WFPC2 camera onboard HST to detect 


( 1999a(b|cl 


detect 102 

addi- 


tional clusters. They derived their physical parameters using in¬ 
tegrated photometry that were compared with simple stellar pop¬ 
ulation (SSP) models, which are ideal star cluster models in the 
sense that they are based on a continuously sampled initial mass 
function (IMF). |Sarajedini et al.| ( |1998| |2000| |2007| ) also used 
WFPC2 and ACS images, but derived the parameters using re¬ 
solved color-magnitude diagrams for the most massive clusters. 
All studies before 2007 have been combined in a merged catalog 
by |Sarajedini & Mancone|p007| . 

More recently, studies in the literature continue to use the 
SSP models to derive the M33 cluster physical parameters (age, 
mass, extinction, and metallicity). However several recent stud¬ 
ies (e.g., |Santos & Frogel|1997[|Deveikis et al.|2008[|Fouesneau| 


& Langon 2010|) have brought attention to the fact that these 
models with continuously sampled IMF, which are unphysical, 
are oversimplified and biased: 


- They are oversimplified because they do not take the natural 
dispersion of star cluster integrated colors into account, the 
so-called stochasticity problem. In reality, the masses of stars 
are stochastically sampled following the stellar IMF, and this 
could be seen as a probability distribution function ( [Santos & 
Frog el|1997| ). Hence, two clusters with same physical param¬ 
eters could host a different number of massive bright stars, 
resulting in very different integrated colors, especially in the 
case of young and low-mass clusters that contain only a few 
massive bright stars. 
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SSP models are biased because they do not even match the 
average of integrated color distributions for star clusters with 
mass log 10 (M/M©) < 4 (see, e.g., [Fouesneau & Langon 


2010t|Popescu & Hanson|2010] [Silva-Villa & Larsen|201lfr 


As a consequence, the derivation of star cluster parameters using 
the SSP models can lead to severe biases, as was demonstrated 
by [Fouesneau & Langon| ( [20101 ) (see their Fig. 3) by use of arti¬ 
ficial tests. 

New models (Fouesneau & Langon|2010| Popescu & Hanson 


|20l0l|da Silva et al.|2012[|de Meulenaer et ai.|20131|2014[|20151 

hereafter Papers I, II, and III) take into account this natural dis¬ 
persion of integrated colors due to the star mass stochastic sam¬ 
pling and allow the cluster parameters to be derived in a more 
realistic way, avoiding the strong biases of the SSP method. 

In this paper, we derive the physical parameters of a merged 
catalog of M33 star clusters, the first time ever using a stochastic 
method on these galaxy star clusters. To this end, we use three 
catalogs of clusters recently published in the literature: the San 
Ro man et al. ( 2010| ) catalog in t he u*g'r'i'z' photometri c system, 


and the |Maf ( |2012| |2Q13| ) and |Fan & de Grijs] ( |2Q14| ) catalogs 
in the standard UBVRI photometric system. We also use near- 
infrared photometry that we derive from deep 2MASS images. 


The structure of the paper is the following. Section[2]presents 
the catalogs of clusters used in this study as well as the 2MASS 
photometry derivation. Section [3] presents the derivation of the 
star cluster physical parameters using our developed method 
which takes stochastic dispersion of star cluster integrated colors 
into account. Section [4] presents the artificial tests used to esti¬ 
mate the reliability of our method for deriving the star cluster 
physical parameters and the global evolutionary fading and dis¬ 
ruption timescales. The derivation of the physical parameters of 
the M3 3 cluster population and the derivation of the fading and 
disruption timescales are performed in Section [5] Conclusions 
are presented in Section [6] 


2. The cluster data 


Recently San Roman et al. (2010) observed 803 M33 star clus¬ 
ters (599 candidates and 204 confirmed using the HST) using the 
3.6 m Canada-France-Hawaii-Telescope (CFHT) and published 
a catalog in the MegaCam camera u*g'r'i'z' photometric sys¬ 
tem. Although their catalog also contained the cluster photome¬ 
try converted to the ugriz photometric system of the Sloan Digi¬ 
tal Sky Survey (SDSS), we considered here the native MegaCam 
photometric system to avoid likely conversion approximation. 

Using archival images of the Local Group Galaxies Survey 
(LGGS, Massey et al.||2006| ) obt ained using the 4 m Kitt Peak 
National Observatory telescope, |Ma| ( |2012| ) built a catalog of 
392 clusters, and |Ma| ( |2013| ) added 234 others , all in the standard 
UBVRI photometric system. |Fan & de Grijs| ( [2014| ) also reana¬ 
lyzed the LGGS photometry to publish a catalog of 708 clus¬ 
ters. In this paper, we adopt for the UBVRI photometric system 
the|Ma|(2012| 2013) photometry when available, and the Fan &| 
de Grijs (2014 ) photometry for the other clusters, correcting the 
Fan & de Grijs (2014) photometric zero-points to the Ma (2012 


2013J) ones in order to have an h omogeneous cata l og. Th e zero- 


point correction coefficients for |Fan & de Grijs | ( |2014| ) minus 
Ma| ( |2012[|2013| ) photometry are AV = -0.099 mag, A (U - V) = 
-0.091 mag, A(B - V) = -0.066 mag, A(V - R) = 0.036 mag, 
and A(V - I) = 0.086 mag, computed for respectively 289, 208, 
208, 250, and 178 clusters with available photometry in common 
between jFan & de Grijs| ( |2014] ) and |Ma| ( |20 1 2| |20 1 3} catalogs. 
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Fig. 1 . The different catalogs of clusters used in this paper: Fan & de 
Grijs ( 201 4]) in open red circles, Ma ( 2012 2013) i n open blue squar es, 


San Roman et al.| ( |2010| ) in green diamonds, !San Roman et al.| (2009 ) ir 


violet crosses. The three dashed ellipses have semi-major axes of 10', 
20', and 30' to the center (marked as a large black plus symbol) and can 
be seen as circles of the same radii projected on the M33 disk. 


These catalogs include all objects published in the catalog of 
star clusters of |Sarajedini & Manconej p007| ), who merged all 
the M33 star cluster catalogs published before 2007. The associ¬ 
ation of these catalogs in this paper results in a merged catalog 
of 910 objects, which is shown in Fig. [T] color-coded to show to 
which original catalog they belong. 

We supplemented optical data with near-infrared data by 
using deep Two Micron All Sky Survey (2MASS) JHK im- 
ageq] with exposure times 6 times longer (2MASS 6X) than 
the standard 2MASS ones. This results in a photometry approx- 
imatively 1.5 mag deeper in 2MASS 6X than photometry de¬ 
rived in standard 2MASS images. The photometry of clusters 
was derived by use of aperture photometry using the standard 
IRAF/digiphot/apphot package with an aperture radius in the 
range 2" to 4", and an aperture correction built using the curves 
of growth of a dozen relatively isolated clusters. By this process, 
we derived the JHK photometry of 758 clusters. To ensure that 
the 2MASS 6X images were correctly calibrated, we also de¬ 
rived aperture photometry of stars, selected in the same region 
where clusters are located, and compared this derived aperture 
photometry to the stellar aperture photometry provided by the 
2MASS Point Source Catalog (2MASS PSC^J which has been 
compiled using the standard 2MASS images. The first row of 
panels in Fig. [ 2 ] presents the comparison of aperture photome¬ 
try of stars derived in this work versus the aperture photometry 
provided by the PSC (black circles) for the JHK passbands. For 
most of stars the agreement is satisfactory. We noted that 109 


1 Kindly made available by T. H. Jarrett, IPAC/Caltech 

2 http://irsa.ipac.caltech.edu 
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Fig. 2. Top row: black circles show the comparison of the aperture photometry of stars derived in this work versus the PSC values, and red 
circles are the comparison of our aperture photometry of the star clusters versus the PSC values for the 109 clusters for which the PSC has 
aperture photometry. Central row: black squares show the photometric uncertainties provided by the PSC for 109 clusters and all circles show the 
photometric uncertainties derived for our aperture photometry for 758 clusters. The red circles are 502 clusters with uncertainties lower than the 
limit of 0.3 mag (dashed line) in all JHK passbands, and the open circles are the clusters for which the photometry does not satisfy this criteria. 
Bottom row: distributions of cluster brightness: white histograms contain the 758 clusters for which we derived photometry and the gray ones 
contain the 502 clusters that satisfy the photometric uncertainty criteria. The thick open histograms represent the 109 clusters contained in the 
PSC, for comparison. In each row, the situation is shown for J (left panels), H (central panels), and K (right panels) passbands. 


clusters are included in the PSC and we could compare the aper¬ 
ture photometry that PSC provides for them with our aperture 
photometry, in red in the figure. The error bars reflect the pho¬ 
tometric uncertainties derived in our aperture photometry. The 
photometric uncertainties of the 758 clusters derived in this work 
are shown in the central row of Fig.|2](all circles), compared with 
the uncertainties of the 109 clusters for which PSC also provides 
aperture photometry (black squares). In this paper we use only 
the clusters with photometric uncertainties lower than 0.3 mag 
in all JHK passbands, indicated in the central row of panels of 


Fig.[2]by dashed lines. This selection results in 502 clusters with 
full JHK photometry. The bottom row in Fig.[2]presents the dis¬ 
tribution of the clusters in each of the JHK passbands of our de¬ 
rived photometry for the whole sample of 758 clusters (in white) 
and for the sample of 502 clusters with uncertainties lower than 
0.3 mag in all JHK passbands (in gray), compared to the distri¬ 
bution of 109 clusters that are included in the PSC (open thick 
histogram). 
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Fig. 3. Color-color diagrams of the cluster sample in optical photometric passbands. The first row of panels shows the color-color diagrams in the 
standard UBVRI photometric system and the second row in the u*g'r'i'z' one. In all panels, the blue circles are clusters that are bright in UV, cyan 
circles are clusters faint in UV, white circles are clusters embedded or close to HII zones, red circles are globular-like clusters, and green stars 
are likely stars rather than clusters. In each panel, the extent and density of the stochastic star cluster model grid is shown as a density surface 
(displayed in logarithmic scale), and the solid line is the SSP model for comparison. In both stochastic and SSP models shown here, the metallicity 
is [M/H] = -0.4. The mass of the stochastic cluster models shown here is fixed to log 10 (M/M o ) = 3.5, a typical mass of the low-mass clusters 
studied in this work, to show the extent of their colors. In each panel, the black arrow shows the Milky Way extinction law direction and the red 
arrow shows the LMC extinction law direction, both computed for A v = 1 mag. 


We also add ultraviolet aperture photometry from GALEX^j 
by using aperture radii of 3"in both far-ultraviolet (FUV) and 
near-ultraviolet (NUV) passbands. This photometry was not 
used to derive the star cluster parameters, but only for the qual¬ 
itative confirmation of results in case of young clusters, as ul¬ 
traviolet magnitudes fade very quickly with age, becoming too 
faint at the distance of M33 after >100 Myr. Also, the very wide 


Point Spread Function (PSF) of the GAFEX telescope makes the 
accurate derivation of UV colors impossible for all clusters, but 
only for the few relatively isolated ones (at least at a distance of 
two aperture radii distance from any other UV-emission in the 
worst cases). 

We also used deep BVRIH a optical images from Subaru 8 m 
telescope and 24 pm image from the Spitzeijj telescope. We cre- 


3 GALEX: https://archive.stsci.edu/prepds/galex_atlas/index.html 
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4 Spitzer: http://ssc.spitzer.caltech.edu/spitzerdataarchives/ 
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Fig. 4. Same as in Fig. [ 5 ] but for GALEX-optical colors (top panels) and optical-2MASS colors (bottom panels). In the top panels, a small part of 
the clusters is shown because most of the clusters are too dim in GALEX passbands, or situated too close to bright neighbor objects. 


ated multi-passband images for each cluster in our catalog that 
were used to visually confirm the results of our method of star 
cluster parameter derivation described in Section [3] 

Figure]?] presents the multi-pas sband photometric data in 
different color-color diagrams in optical cases (UBVRI and 
u*g'r'i'z'). Figure]?] shows the clusters in GALEX photometry 
(top panels) only for objects undisturbed by close neighboring 
UV emission, and in deep 2MASS 6X photometry (bottom pan¬ 
els). The clusters are shown with SSP models (solid lines) and 
also with a grid of artificial star cluster models, which take the 
stochastic dispersion of their colors into account, as described in 
Section]?] Both SSP model and stochastic model grid are shown 
with the same metallicity, [M/H] = -0.4. 

Although GALEX photometry is inaccurate for most of the 
910 objects because of the presence of possible UV emitting 


neighboring objects, we used the FUV photometry as a crite¬ 
rion for the qualitative evaluation of the UV emission strength. 
Objects are said to be bright in UV when their aperture photom¬ 
etry within an aperture radius of 3" is brighter than 20 mag, and 
faint in UV when it is fainter than this limit. 


In Figs.]?] and [4] the clusters are color-coded according to 
their types: bright in UV (blue), faint in UV (cyan), embedded 
or close to HII zones (white), and globular-like clusters (red). 
Clusters are classified as globular-like depending on their visual 
confirmation usi ng HST images ([Sarajedini et al.||1998[ |Chan- 

dar et al.||1999a[ |San Roman et al.||2009| ). In addition, the deep 

Subaru images were used to reject 95 highly probable stars from 
our star cluster sample. These stellar objects are marked as green 
star symbols in the color-color diagrams in Figs.]?]and]?] 
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After the rejection of these 95 stars from the 910 objects in 
the sample, as well as a few clusters with incomplete photometry, 
747 clusters remain to be studied. 


3. Method of derivation of cluster parameters 


Following Fouesn eau & Lan gon (2010|); [Fouesneau et al.| ( |2014] ), 
and Papers I, II, and III, the derivation of the physical parameters 
(age, mass, extinction, and metallicit)^) of a given observed star 
cluster is based on a comparison of its integrated broad-band 
photometry to a four-dimensional grid (for the age, mass, ex¬ 
tinction, and metallicity) of star cluster models. Each node of the 
grid contains 1 000 star cluster models. Each star cluster model is 
built by randomly sampling the stellar mass according to the IMF 
( Kroupa|2001| ) following the method desc ribed in[D eveikis et al 
(2008]) (see also |Santos & Frogel|[l997| |Cervino et aLpOO^f . 
The luminosities of clusters were der ived using the PADOVA 
isochrone^] from Marigo et al. (2008 ) with the addition of the 
Thermally Pulsing Asymptotic Giant Branch (TP-AGB) phase 
from |Girardi et al.| ( [2010 ). The grid was built according to the 
following nodes: from log 10 (£/yr) = 6.6 to 10.1 in steps of 0.05, 
from log 10 (M/M o ) = 2 to 7 in steps of 0.05, and for 13 metal- 
licities: from [M/H] = +0.2 to -2.2 in steps of 0.2. This results 
in a grid of 71 values of age, 101 values of mass, with 1 000 
models per node, hence ~7 x 10 5 6 models for each metallicity. 
To limit the number of models that need to be stored in com¬ 
puter memory, the extinction was computed when the observed 
cluster was compared with the grid of models. It ranges from 
E(B - V) = 0 to 1 in steps of 0.01, therefore 101 values for the 
extinction. We used the [Gordon et al.| ( [2003] ) extinction law de¬ 
rived for the Large Magellanic Cloud (LMC), as the M33 galaxy 
is believed to have a similar metallic content and hence may fol¬ 
low a similar extinction law. 

We evaluated the likelihood of each node of the grid to rep¬ 
resent the magnitudes of a given observed cluster. Within each 
node, we first computed the likelihood of each of the 1 000 star 
cluster models by 


Gnodel — /- 

/=! V2 ncTf 


exp 


(magy obs inagy-model) 


2 °"/ 


( 1 ) 


where / stands for one particular filter, mag^ for the observed 
and model magnitudes in that filter, and F for the total number of 
filters. For example, F = 5 for the UBVRI photometric system 
we use in this study. Then the likelihood of the node of age t , 
mass M, extinction E(B - V), and metallicity [M/H] is the sum 
of the likelihoods of its models, 


N 

L a ode (t, M, E(B - V), [M/H]) = 2 L model ,„ , (2) 

n= 1 


where N = 1 000, the total number of models contained in 
the node. The procedure is repeated for each node of the four¬ 
dimensional grid, and the observed star cluster is then classified 
with the parameters of the node, which maximizes the quantity 
L node . We note that this procedure could also be applied by using 
colors (e.g. U - B, u* - g', or other passband combinations) in 
place of individual magnitudes as the variable mag^ of Eq. [T] 

5 Hereafter we refer to extinction and metallicity as E(B - V) and 
[M/H]. 

6 PADOVA isochrones from “CMD 2.6”: http://stev.oapd.inaf.it/cmd 
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4. Artificial tests 

The ability of the method to derive star cluster parameters has 
already been evaluated in Papers I, II, and III. Here we are in¬ 
terested in seeing for which conditions the method would be 
sensitive enough to detect a change of slope in the number of 
clusters per age bin distribution (hereafter referred as differen¬ 
tial age distribution ) of the cluster sample such as shown by the 
solid line in Fig.[5]i Indeed, a two-slope profile in the differential 
age distribution could be interpreted as a decrease in the number 
of clusters due to an evolutionary fading of the cluster magni¬ 
tudes (first slope), and a decrease of the number of clusters due 
to their disruption (second slope), as is discussed in greater de¬ 
tails in Sect. [5] Here the objective is to model an artificial star 
cluster population with such a two-slope profile in the true dif¬ 
ferential age distribution and to see whether the derived differen¬ 
tial age distribution reproduces this profile depending on which 
photometric system we use. 

We generated a sample of 10000 artificial star clusters. The 
differential age distribution of the artificial clusters was chosen 
to mimic a two-slope profile similar to that described in |Van- 
|sevicius et al.| ( |20 09 ) for M31 star clusters (see their Fig. 5a, 
reproduced in our Fig.[5ji in solid line). For simplicity, the mass 
of the clusters was fixed to log 10 (M/M o ) = 4, a typical value 
for the clusters observed in M33. The extinction was randomly 
generated uniformly in the range E(B -V) = 0 to 1. 

For the cluster metallicities, we use a very simple age- 
metallicity relation: for the youngest clusters the metallicity is 
supersolar, [M/H] = 0.2, and for the oldest clusters the metal¬ 
licity is very low, [M/H] = -1.8. The age-metallicity relation is 
linear between these values in the age (log 10 (r/yr))-metallicity 
([M/H]) space. The metallicity of clusters is generated with a 
Gaussian dispersion of 0.4 dex standard deviation around this 
age-metallicity relation. 

The first test, presented in Fig.[5j was performed using a 
cluster model grid with the metallicity fixed to [M/H] = -0.4 
to show the possible biases that occur when metallicity is 
not a free parameter. The test was performed for three pho¬ 
tometric system combinations: optical (UBVRI), optical with 
near-infrared (UBVRI + JHK) and ultraviolet with optical 
(GALEX+ UBVRI). The UBVRI system was used as a refer¬ 
ence for the mass derivation, so it was used in magnitudes in 
the Eq.[l] For the other passbands used we have used colors in¬ 
stead, FUV-NUV for the GALEX passbands and J - H, J - K , 
and H - K for 2MASS. This allowed us to combine different 
catalogs of clusters built with slightly different aperture sizes: 
we used the magnitudes for one catalog, and the colors for the 
others. However it is important that for each cluster at least one 
magnitude should be given, not just colors, so that the mass of 
the clusters could be estimated reliably by the method. 

We added photometric uncertainties as a Gaussian noise with 
standard deviations of 0.05 mag for each UBVRI photometric 
passbands, 0.1 mag for JHK , and 0.15 mag for GALEX FUV 
and NUV. 

In Fig. [5] we see how the derived cluster age (panel a) and 
mass (panel b) are distributed versus the true values (indicated by 
dashed lines). Fig. [5]: and Fig.[5]l concentrate on the age deriva¬ 
tion, and show that the peak in the true age distribution (indicated 
by thin lines in panels c and d) has already been found using op¬ 
tical data only. The true peak in age, situated at log 10 (f/yr) ~ 8.3 
in Fig. [5]: (solid line histogram) and corresponding to a change 
of slope at log 10 (f/yr) = 8.5 in Fig.BH (solid line) is correctly de¬ 
rived as maximum in Fig. [5]: (dashed line histogram) and change 
of slope in Fig.[5]l (dashed line). However, the apparent good 
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Fig. 5. Artificial tests for 10000 clusters with true age and metallicity following a defined age-metallicity relation (see text), and with a true 
mass fixed to log 10 (M/M o ) = 4. The first column of panels displays the age derived vs the true age, the second column shows the derived mass 
distribution, the third column shows the true age distribution (solid line histogram) and derived age distribution (dashed line histogram), and the 
last column shows the differential true age distribution (solid line), which is composed of a two-slope profile, and the derived differential age 
distribution (dashed line). The first row of panels shows results obtained using the UBVRI photometric system, the second row shows the results 
obtained using the UBVRI + JHK system and the last row shows the results obtained with the GALEX + UBVRI system. Here the metallicity of 
the model grid is fixed to [M/H] = -0.4. 


match between true age and derived age distributions in Fig. |5j: 
and Fig.[5}l can be rather misleading when using UBVRI pho¬ 
tometry alone. Indeed the direct comparison of the individual 
clusters’ true and derived age in Fig. [5^ shows that the agreement 
is far from evident, especially at old ages, where the true metal¬ 
licity of artificial clusters and the fixed value of the model grid 
([M/H] = -0.4) deviate most. As a natural consequence, most 
of the clusters with true age above log 10 (f/yr) = 9.5 have age 
underestimated. Also, the presence of the natural age-extinction 
degeneracy in the optical UBVRI case, already discussed in pa¬ 
pers I and II, produces the streaks developing perpendicularly to 
the left and to the right of the diagonal identity dashed line in 
Fig. [5^. The situation is less extreme, but still strongly affected 
by these degeneracies, when we add JHK passbands to UBVRI 
ones (second row of panels). When we use GALEX with UBVRI 
passbands (third row of panels), the age-extinction degeneracy 
disappears, but the deviation from the identity line occurs be¬ 
cause of the strong sensitivity of ultraviolet to the metallicity. 
[Bianchi[f2011| indeed shows by use of integrated spectra of sim¬ 
ple stellar population models that it is in the ultraviolet spectral 
region that the spectra are most affected by a change in metallic¬ 
ity. 


We performed a second test, fixing the metallicity to 
[M/H] = -0.4 for all clusters, and then, only for clusters which 
have derived age larger than log 10 (£/yr) = 9, we re-derive a solu¬ 
tion leaving the metallicity free to vary in the range [+0.2, -2.2]. 
Indeed one notices in the first test that the situation was most 
complicated for clusters with true age above 1 Gyr. The results, 
shown in Fig. [6] still suffer from strong age-extinction degeneracy 
in the case of UBVRI passbands only (first row of panels). The 
inclusion of near-infrared photometry improves much the deriva¬ 
tion as the streaks developing perpendicularly to the identity line 
in Fig. [6^ are strongly reduced. In this case, the match between 
the true and derived age distributions in Fig. [6^ and Fig.[6ji is 
much more secure for all age ranges. When using GALEX with 
UBVRI (third row of panels), a gap is visible in Fig. [6} due to 
the strong sensitivity of the ultraviolet flux to metallicity. In this 
case, strong biases are still present in the age distribution (Fig.[6jc 
and Fig.[6}). 

A third test was performed in which metallicity of the model 
grid was left free in the whole age range, and the results are pre¬ 
sented in Fig. [ 7 ] Here we see that the use of optical passbands 
only (first row of panels) or even optical with near-infrared (sec¬ 
ond row) can lead to strong biases as these photometric systems 
are not sensitive enough to discriminate between models of dif- 
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Fig. 6. Same as in Fig.[5j but here, for clusters which have derived ages larger than log 10 (f/yr) = 9 when using fixed [M/H] = -0.4 metallicity, we 
re-derive a solution leaving the metallicity free to vary in the range [+0.2, -2.2]. 


ferent metallicities (see also Papers II and III for the sensitivity 
of the derived parameters on the metallicity, as well as for the 
derivation of the metallicity parameter). As a consequence, age 
distributions (panels c and d for UBVRI case, panels g and h 
for UBVRI + JHK case) are strongly affected. Only ultraviolet 
associated with optical data passbands are able to break the age- 
extinction degeneracies when metallicity is left free, as shown in 
last row of panels. As a consequence, the derivation of the cor¬ 
rect two-slope profile in the differential age distribution is best 
done using GALEX + UBVRI when metallicity is left free, see 

Fig.0)- 


5. Derived physical parameters of clusters 

We applied the method of derivation of physical parameters 
to the sample of 747 M33 clusters using the optical UBVRI 
and near-infrared JHK passbands. We first fixed the metallic¬ 
ity to [M/H] = -0.4 for all clusters, and then, only for clus¬ 
ters that have derived ages larger than log 10 (f/yr) = 9, we re¬ 
derive a solution leaving the metallicity free to vary in the range 
[+0.2, -2.2], as was shown to be the best choice for this pass- 
band combination in the previous section. 

As was done for the artificial tests, we used the UBVRI sys¬ 
tem as a reference for the mass derivation, and so it was used in 
magnitudes in Eq.[l] For the other passbands used, u*g'r'i'z' and 
JHK , we used colors instead to avoid problems of different aper¬ 
tures in the different catalogs used. The colors used were u* -g\ 
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g' - r', g' - i', r' - i', and i' - z! for the CFHT passbands, and 
J - H, J - K, and H - K for the 2MASS passbands. 

We used the extinction law of [Gordon et al.| ( [2003| ) derived 
for the EMC, assuming that for a similar metallic content the 
M3 3 galaxy would have a similar extinction law. The minimum 
extinction of clusters was set to E(B -V) = 0.04 mag, the value 
of the foreground extinction in the direction of M3 3 estimated 
from the Schle geTet al.| ( |1998] ) extinction maps. 

The results obtained here are compared in Fig.[8^,b,e,f for the 


age and mass with the ones of 160 clusters of San Roman et af 


(2009), obtained by isochrone fitting on HS T-resolved color - 
magnitude diagrams. In the case of the |San Roman et aL] ([2009), 
cluster ages are enclosed in a much narrower range, mainly be¬ 
tween 50 Myr and 1 Gyr. Although our age distribution is wider 
than in their case, a satisfactory agreement is found between both 
sets of results, as well as for the mass parameter. 

Globular-like clusters (red circles, visually confirmed as 
globular clusters on HST images) are found to be very old in 
our case, and more massive. |San Roman et al.| ( |2009| ) noted that 
the lack of clusters with ages older than ~ 1 Gyr in their cata¬ 
log is linked to the fact that their resolved color-magnitude di¬ 
agrams are generally not deep enough to detect the main se¬ 
quence turn-off of clusters older than this age. In Fig.[8^,b,e,f, 
the two oldest and most massive clusters (according to our de¬ 
rived age and mass) are globular clusters known as CBF85-U137 
and MKKSS12. For CBF85-U1 37, we find log 10 (i/yr) = 10.05 
and log 10 (Af/Af o ) = 5.25 while San Roman et al. (|2009|) give 
logio^/yr) = 8.8 and log 10 (M/M o ) = 4.4. Chandar et al.| ( |2Q02 ) 
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Fig. 7. Same as in Fig. [5] but here the metallicity of the model grid is left free over the whole age range. 


give log 10 (f/yr) = 10.08 (12Gyr in their table 5) using spec¬ 
troscopic line index models of | Worthey] ( | 1994| . For MKKSS12, 
we find log 10 (f/yr) = 10.0 and log 10 (M/M o ) = 5.65 while |San 
Roman et~ak| (2009 ) give lo g 10 (f/yr) = 8.6 and log 10 (M/M o ) 


4.54. [Chandar et ah ( 2002|) giv e log 10 (f/yr) = 9.4 using SSP 
models, and |Ma et al.| (|2004) give log 10 (?/yr) = 9.63 and 
log 10 (M/M o ) = 5.47 using their BATC spectrophotometric sys¬ 
tem composed of 13 narrow passbands, also compared to SSP 
models. 

For young clusters, the age given by San Ro man et al.| (2009) 
is often older than in our values. In Fig.[8^,e we note that the two 
white circles, which are clusters still within HII zones, and so 
very young, are also seen as older in |San Roman et al.| ([2009 ) 
than in our study. Also, many clusters that are bright in UV 
(blue circles) are also older in |San Roman et al.| ( |2009| ) than in 
our case. However, UV brightness fades very quickly, becoming 
faint at the distance of M3 3 after >100Myr, making it unlikely 
that the age of these clusters is older. 

Kg© ,d,g,h compares the results found for clusters com¬ 
mon with the |Sarajedini & Mancone|p007| ) merged catalog (la¬ 
beled “SM10” in the figure because revised in 2010); the agree¬ 
ment for both age and mass parameters is not as good. One has 
to keep in mind that the [Sarajedini & Manconej ( |2007| ) cata¬ 
log contains results from very different studies. Many of the¬ 
ses results have been obtained using SSP models on integrated 
unresolved ground-based photometry using optical colors only, 
a technique that has been shown to be strongly biased by the 
stochasticity problem and the presence of strong degeneracies 
(see, e.g., [Fouesneau & Langon|2010[ Papers I and II). Foues- 


neau & Langon ( |2010| ) created a sample of stochastically gener¬ 
ated artificial star clusters and tried to derive their age using SSP 
models. Their results, shown in the bottom left panel of their 
Fig. 3, are similar to the values given in our Fig.[8]:. The concen¬ 
trations of solution at ages log 10 (f/yr) ~ 7 and log 10 (r/yr) ~ 9 
for |Sarajedini & Mancone| ( [2007| ) seem to be artifacts due to the 
S SP method, as is the case for the SSP derived ages in Fig. 3 
of Fouesneau & Langon (2 0T0| ). |Popescu et al.| ( |2012] ) show the 
same features, this time for real clusters. They derived the age 
of LMC star clusters using the stochastic method and compared 
these values to the ages derived by |Hunter et al.| ( [2003]) using the 
SSP method. The comparison, shown in Fig. 8 of |Popescu et al. 
( 2012[ ) (the order of the x- and y-axes is flipped co mpared to 
our figure) shows similar features to the Fouesneau & Langon 
( |2010 ) study and to this work: the deviations from the identity 
line are attributed to artifacts of the SSP method of parameter 
derivation. 

Figure[9] presents the results for the star cluster sample stud¬ 
ied in the paper. As expected, the mass versus age distribution 
in Fig.[9^ shows a different typical age for the different classes 
of clusters. The youngest ones are the embedded or close to HII 
zones (white circles), then come the clusters bright in the UV 
(blue circles), then the faint in the UV (cyan circles), and finally 
the globular-like clusters (red circles). The solid line shows the 
50% completeness line evaluated in |Fan & de Grijs| ( |2014| ) for 
their cluster sample. As our merged cluster sample also contains 
HST detected objects from [San Roman et al. ( 2009]), some clus¬ 
ters may be found well below this limit. 
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For M33, 
giant stars, w 
on E(B - V ) 


U et al.| ( |2009| ) derived the extinction of 22 super- 


lich resulted in an extinction distribution centered 
0.1 mag . |U et al.| ( [2009) ) also used the data of 
[Rosolowsky & Simon] ( 2008| T to derive E(B - V) values for 58 
HII regions, and show in their Fig. 9 that extinction can be ex¬ 
pected to be E(B - V ) < 0.3 mag for those regions (except for 
three objects), with an average E(B - V) ~ 0.11 mag. 

The extinction of all clusters depending on their deprojected 
galactocentric distance (assuming that all clusters are in the disk, 
which is incorrect for globular-like clusters) is shown in Fig.[9]f. 
The global median extinction is 0.16 mag. The extinction that we 
found is generally higher for young clusters than for old ones, as 
shown in Fig. (9]}. Indeed, the majority of clusters still embedded 
or close to HII zones (white circles) are found to be more ex- 
tincted with a median of 0.34 mag. Clusters bright in UV (blue 
circles) have a median extinction of 0.17 mag while clusters faint 
in UV (cyan circles) have a lower median extinction of 0.14 mag. 
Globular-like clusters (red circles) have the smallest median ex¬ 
tinction with 0.09 mag. 

Figure[l0| describes the age and mass distributions (panels 
a and b) as well as the differential age and mass distributions 
(panels c and d). We see that the differential age distribution is 
com posed of a two-slope profile. [Boutloukos & Lamers ( 2003) 
and |Lamers et al.| ( |2005| ) interpreted the first slope as a natural 
magnitude fading as a result of stellar evolution, and the sec¬ 
ond slope as being due to cluster disruption mechanisms such 
as the galaxy tidal field effect or encounters with giant molec¬ 
ular clouds. Hence we see here that the cluster sample is dom¬ 
inated by the magnitude fading until log 10 (f/yr) ~ 8.5 and that 
after the cluster disruption phase takes over. This typical life¬ 
time scale is comparable to that derived for the star cluster pop- 
ulation in th e southwest field of the M31 galaxy by Vansevicius 


et al.|(2009]) usin g a star cluster sample photometry from Nar- 
butis et al.|(|2008|). 


Following Gieles | < [2Q09|) , we co mpare the cluster differential 
mass distribution to the |Schechter| \\916) distribution function, 
defined as 


dN/dM = A x M~P x exp(-M/M*) 


(3) 


where A constant scales with the cluster formation rate, is the 
power-law index of the mass function and M* stands for the char¬ 
acteristic mass after which the exponent term decreases strongly. 

As was found for most spiral galaxies (see, e.g., |Larsen] 2009; 
Portegies Zwart et al. 2010| ), the derived mass function of the star 
cluster sample follows the Eq. [3] distribution function with J3 = 2 
and M* = 2 x 10 5 M o . We adapted the scaling constant A to scale 
it to the cluster mass distribution, shown in Fig. ITOli. 


6. Conclusions 

We studied the star cluster system of the M3 3 galaxy, using the 
most recent optical broad-band photometry catalogs, and sup¬ 
plemented near-infrared measurements using deep 2MASS im¬ 
ages. As most of clusters are partially resolved or unresolved, we 
used a method of star cluster parameter derivation which takes 
into account the natural dispersion of the integrated colors due 
to the stochastic sampling of stars in the clusters. We present the 
derivation of the age, mass, and extinction of the clusters for a 
metallicity fixed to [M/H] = -0.4 (LMC-like), and when the age 
derived is larger than 1 Gyr, then a new solution is derived using 
free metallicity in the range [+0.2, -2.2]. 

We ensured, by use of artificial clusters, that the star clus¬ 
ter physical parameter derivation method can correctly derive a 

Article number, page 10 of 


no. Article_IV_M33_final 



b) 


20 

80 ■ M 

40 +1 
J n-TI-Ln 111 _ 

2 3 4 5 6 

log lo (M/M 0 ) 



Fig. 10. Top row: age (panel a) and mass (panel b) distributions derived 
for the M33 star cluster sample. Bottom row: differential age (panel c) 
and mass (panel d) distributions. In panel c, the solid line and the dashed 
line are respectively the cluster evolutionary fading rate and the cluster 
disruption rate, both taken from the case of M31 ( [Vansevicius et al.| 
|2009| and shifted vertically here, for comparison. The dark histogram 
in panel b and the black circles in panel d represent a subsample of 
clusters with ages between 100 Myr and 3 Gyr. The dashed line in panel 
d represents the cluster mass function that follows a |Schechter||l976| 
function with J3 - 2 and M* - 2 X 10 5 M o . 


given two-slope profile in the differential age distribution, test¬ 
ing it for different photometric systems: optical alone ( UBVRI ), 
optical with near-infrared (UBVRI + JHK ), and ultraviolet with 
optical (GALEX+ UBVRI). We showed that the optical with 
near-infrared case is fit for the correct derivation of the two-slope 
profile, and we used it for the M33 star cluster system. 

A two-slope profile of differential age distribution shows that 
the typical lifetime before disruption of star clusters in the M3 3 
star cluster system is found to be ~ 300 Myr, comparable to 
what is found for M31 star clusters. We show that the differ¬ 
ential mass distribution of clusters is consistent with a lSchechterl 
( 1976| function with a power-law index J3 = 2 and a characteris¬ 
tic mass M* = 2 x 10 5 M o . 
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Fig. 8. Comparison of the age and mass derived for clusters common with San Roman et al. < |2009| (panels a,b,e,f), and for clusters common 
with Sarajedini & Mancone ( 2007) (panels c,d,g,h). The first row shows the comparison of the parameters of these studies vs this work, and the 
second row shows the difference of the parameters derived in their studies and this work vs the parameters derived in this work. The color-coding 
of clusters is the same as defined in Fig.[3] 




log 10 (f/yr) E(B—V) deprojected distance [arcmin] 


Fig. 9. The first row of panels shows the mass vs age (panel a), extinction vs age (panel b), and extinction vs mass (panel c) derived for the 747 
clusters studied in this paper. The solid line in panel a shows the photometric 50% completeness limit estimated in |Fan & de Grijs| < |2014|) for 
their optical ground-based cluster catalog. The second row shows the metallicity vs the age (panel d), the extinction histogram (panel e) and the 
extinction vs the deprojected galactocentric distance (panel f) of the cluster population. As emphasized in the text, the metallicity of the clusters, 
shown in panel d, has been fixed to [M/H] = -0.4 when the derived ages using that metallicity is lower than 1 Gyr, and for clusters with older 
derived ages, a new solution is derived with free metallicity in the range [+0.2, -2.2]. The color-coding of clusters is the same as defined in Fig.[3] 
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